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ABSTRACT 






Image theory is an ideal method for calculating the transmission loss in a shallow 






water (wedge shaped ocean) environment. It can be used in cross-slope, at all 






frequencies and in transitional cut off regions that are out of bounds to normal mode 






theories. 
This thesis had three objectives: 1) convert the existing image theory models called 
URTEXT and WEDGE into a high level scripting language called MATLAB™ by 







Math Works, 2) linearize the existing quadruplet expansion program to increase speed, 






and 3) to incorporate the Arctan approximation of the Rayleigh reflection coefficient 






into the quadruplet expansion for the fast bottom case. 
Objective 1 was completed with accurate results. Objective 2 was completed with 






a factor of 8 increase in speed. Objective 3 incorporated the Arctan approximation of 







the reflection coefficient for a fast bottom into the quadruplet expansion, but due to the 






inaccuracy of the reflection coefficient after the second quadruplet, the results were not 






favorable. It was also discovered that even with the Rayleigh reflection coefficient, the 
first order approximations made in developing the quadruplet expansion equation 









(Equation 6-27) are not accurate enough for the fast bottom case. 
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I. INTRODUCTION 


The current emphasis of undersea warfare is on the shallow coastal regions in littoral 
waters. The sound propagation in these regions is influenced by bottom and surface 
interaction and not by the convergence zones and ducting found in deep water. The 
variables in shallow water are bottom type, bottom density, slope, sea state, surface winds, 
sea surface temperature, etc., all of which can change quickly or be unavailable in the area 
of interest. 

There is no closed form analytic solution to the problem of a penetrable, sloping 
bottom as found on the continental shelf region or shoreline areas. There are computer 
models using various approximation methods ranging from purely numerical (parabolic 
equation) to approximately analytical (adiabatic mode theory). Transmission loss 
modeling in a wedge shaped ocean falls into three major areas: 1) parabolic approximation, 
2) adiabatic normal mode theory, and 3) image theory. 


A. PARABOLIC EQUATION 
The parabolic equation (PE) is a range-dependent, underwater acoustic propagation 
model. In it the Helmholz equation is replaced by a one-way parabolic equation which 


generates an acoustic field as an initial value problem. The PE approach is limited by 


excessive computer running time in shallow water environments, higher frequencies, and 
horizontal rays < 40°. 

The first application of PE to underwater acoustics was by Tappert in 1977 with a 
restriction on the maximum angular aperture of +20°at the source [Ref. 1]. This restriction 
was relaxed to +40° with a higher order parabolic equation by Collins {Ref. 2]. 
Computerized calculations by Jensen and Kuperman [Ref. 3] based on a split-step Fast 





Fourier Transform (FFT) solution of the parabolic equation provided detailed information 
about the beams projected into the bottom, showing modal cutoffs. The results compared 
favorably to model tank experiments by Coppens and Sanders [Ref. 4]. Since the PE is 
cylindrically symmetric, sloping bottoms cause a problem. Lee and McDaniel [Ref. 5] 
sectioned the. wedge into a series of range independent regions and then applied normal PE 
methods. Another approach to this problem is by Collins [Ref. 6] with the rotated PE. 
The rotated PE steps parallel to the ocean bottom and has two normal derivatives 
preserving pressure and particle velocity normal components thus handling sloping 
interfaces properly. Fawcett [Ref. 7] developed a three-dimensional computer mode! 
using the wide-angle PE approximation by Thomson [Ref. 8] with FFT algorithm as an 
azimuthal operator. Most recently, Collins [Ref. 9] has developed the energy conserving 
PE for elastic media to improve accuracy in range-dependent elastic media. The energy- 
conserving elastic PE is a generalization of the energy conserving acoustic PE [Ref. 10]. It 
involves approximating a range-dependent waveguide into a sequence of range- 
independent regions with a linear approximation of compressional energy flux between the 


vertical interfaces. 


B. NORMAL MODE THEORY 
Normal mode theory provides an exact solution to the wave equation. The distinct 
advantage of normal modes is that once the set of eigenfunctions have been determined, the 
range and depth dependence of the transmission loss can be calculated directly [Ru*. 11]. 
Normal mode theory is a range-independent approach but for a wedge shaped ocean, a 
range-dependent approach is needed due to the bottom interaction (in the case of a lossy 
bottom) and bottom reflection angle. Pierce [Ref. 12] used adiabatic separation of depth 


and range in the wave equation to get an approximation. Graves, Nag], Uberall, and Zarur 
[Ref. 13] applied this method to the wedge shaped ocean using rigid bottom and isospeed 





water. This method was good only for small slopes. Buckingham [Ref. 14] developed a 
solution for the penetrable bottom using an “effective wedge with pressure release bottom 
below the actual penetrable bottom. By using a coordinate transformation, range- 
dependent modes are matched to range independent eigenfunctions in the "effective" wedge 
to account for amplitude and phase shifts. Then a sum of range dependent normal modes 
are applied to the new bottom. 

A problem with adiabatic normal mode approximation is it does not explain the 
transition to evanescent modes at cutoff. Pierce [Ref. 15] combined PE with adiabatic 
normal mode theory to develop a critical depth function. Known as the augmented 
adiabatic mode theory, it is a sum of modal terms with each mode eventually encountering 
a critical depth at a critical range. The results compare well with the results of Jensen and 
Kuperman [Ref. 3]. 

Another method is to combine local modes with ray acoustics. Amold and Felsen 
[Ref. 16] developed intrinsic wave functions that are determined by waveguide geometry. 
Modes are mutually decoupled thus propagate independent of source geometry (physical 
shape) when superpositioned. 


C. IMAGE THEORY 

Image theory in acoustics is similar to image theory in optics. A source radiating in 
the water has a virtual image reflected equidistant above the water 180° out of phase. In an 
isospeed, wedge shaped ocean numerous images are produced in kaleidoscope fashion. 

In 1966, Macpherson and Dainteth [Ref. 17] proposed a phase incoherent model for 
upslope propagation . Kawamura and Ioannou [Ref. 18] predicted the pressure amplitude 
and phase of the sound field along the bottom of a wedge shaped fluid layer overlapping a 
fast fluid bottom. Coppens, Humphries, and Sanders [Ref. 19] developed a phase coherent 
model to calculate the pressure amplitude along the bottom of a wedge of water overlying a 





fast fluid bottom in the upslope direction using the incident plane wave reflection 







coefficient stated by Rayleigh. The pressure distribution at the water-bottom interface can 





be used as a source to project beams into the bottom to calculate the propagation through 






the medium. Brekhovskikh and Lysenov [Ref. 20] stated that for very distant source and 


image combinations or the presence of many reflections allow image theory to be applied 







with Rayleigh reflection coefficients with insignificant error. Baek [Ref. 21] then 






developed a computer model of pressure throughout the water column over a fast botto: 






This model was limited to a fast bottom due to the critical grazing angle of the sound to 






bottom. LeSesne [Ref. 22] used the same model as a basis for a three-dimensional model 
which produced a pressure field in the cross slope case which was validated against 
experimental results from a model tank by Kosnick [Ref. 23]. Kaswandi (Ref. 24] 







developed a model for a pressure field with a slow bottom in downslope. Nassopoulos 






[Ref. 25] took all previous models developed at the Naval Postgraduate School and 






combined them into a computer program called URTEXT. He also implemented the use 
of the doublet expansion of the acoustic field in a wedge shaped ocean. Livingood (Ref. 





26] then expanded the image doublets into the cross-slope case. Joyce [Ref. 27] then 







developed the quadruplet expansion of the acoustic pressure field in the wedge shaped 






ocean. The quadruplet expansion is a much faster computationally but the approximations 






in the reflection coefficient and other terms in the quadruplet expansion equation produces 


large errors at large wedge (B>3°) and source angles due to the angular differences between 





the images in the upper and lower doublets. 





Il. IMAGE THEORY 


This section is a brief review of the ideas and basic equations of image theory. The 


following assumptions are made: 
1. Sound velocities and densities are constant in both the wedge and penetrable bottom. 
2. Air-water boundary is a pressure release surface i.e. reflections are 180° out of phase. 
3. Slope is constant and all boundaries are planar. 


When sound is transmitted within a waveguide, there are a number of interactions of the 
acoustic waves with the boundaries. In image theory, each ray path is replaced by an image 
of the source. The distance of the source from the receiver is the total length of the acoustic 
ray. In the case of a sloping bottom, the image is placed equidistant from the reflecting 
boundary perpendicular from the source. See Figure 2-1. 





Figure 2-1 -Wedge Geometry 
For a given bottom slope angle B, the number of images N in each of the upper and 


lower half spaces is given by 


B 
Each image is numbered from the source. The source is image number | for the upper 


N= inf (2-1) 


half space. The first image below the bottom is image number 1 for lower half space. This is 
continued until the opposite surface is reached. See Figure 2-2. 





Figure 2-2 - Image Structure for Wedge Shaped Duct 
The grazing angle 6, of the nth image above or below the bottom is calculated from 


6.=nB-y forn odd (2-2) 


0, =(n-1)B+y for n even (2-3) 


where ¥ is the angle of the source from the surface. 


The range r,, from each image to the receiver is then calculated using 


r., =r +re + ye —2rr,cos(0, —B +5) 


for the upper series of images and 


r, =r? +r2 + y2 -2r,r, cos(0, +B - 5) (2-5) 
for the lower series of images. In these equations, r; is the distance from the wedge apex to 
the source, 2 is the distance from the apex to the receiver, yo is the cross-slope or shoreline 
distance between the source and receiver, and @,-B+6 is the angle from the nth source to the 


receiver (for the upper image). See Figure 2-3. 


Figure 2-3 - Source-Receiver Geometry (Upper Image) 
The sound will propagate from the source to the receiver through combinations of three 
different path types: 1) direct, 2) surface reflected, and 3) bottom reflected. Figure 2-4 
illustrates a bottom-surface-bottom path. The upper drawing illustrates the real path of sound. 





Receiver 





Figure 2-4 - Comparison of Real and Image Theory Sound Propagation Paths 
The surface is a pressure release boundary so the incident angle is equal to the reflected angle. 
The bottom has refractive properties due its composition so the incident angle is not 
necessarily equal to the transmitted angle but is dependent on the p and c ratio of the water 
and bottom. With each interaction with the bottom, a new angle is produced which is related 
to the first incident angle. This is discussed later with the reflection coefficient. The lower 
drawing illustrates the image theory equivalent. The path from the third image to the receiver 
is a straight line. 





The reflection coefficient is a product of the interaction between the boundaries of the 
fluid layers. The surface is a pressure release boundary so each interaction produces a 
reflection coefficient of -1. Each interaction with the bottom requires calculation of the 
reflection coefficient. This coefficient is a function of the speed of sound in each medium, the 
density of each medium, and the grazing angle of the ray equivalent on the bottom. The sine 
of the grazing angle O,,. is 


eg ah sin(0, — 2mB) + r, sin{(2m — 1)B + 5] 
am+ r 


ar 


for the upper images and 


cy ees sin(6, — 2mB) + r, sin{(2m + 1)B - 5) 
di r 


a- 


for the lower images. See (Figure 2-5). 
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lh 


APEX 


Figure 2-5 - Incident Angle 6,,, Calculation (Upper Half Space) 
The reflection coefficient for each interaction R,, then is given by 





sin 0, - (Ve +a’ +a)" + j5(Vor +a? - a)" 


Fe al ee ee ee 
Pz . Ji 2 2 % Py 1 2 2 A 
7 sin 8,,. + (Vo +a +a) ~ ize (vo +a -a) 


where Pp}, P2 are densities of the water and bottom respectively, j = V-1, and a and b are 


(2-8) 


given by 
2 
a= (2) — cos” On (2-9) 
7) 
him =) a (2-10) 
c,) k, 


For the above equations, c;, c2 are the speeds of sound of the water and bottom, 
respectively, and o/kz is the bottom loss coefficient. 


The pressure for each half space is given by summing the contribution for the upper and 


lower images respectively 
N M 
B= J exp(— jhe TT Rae (2-11) 
azl "nt m=) 
where [] R=1 when n=1,2 in (2-11) and 
ol , wr) 
P, = Y—exp(-Jer,-X-1)"" TR... (2-12) 
n=1 " a- m=0 


The number of interactions m with the bottom is related to the image number, so the limit M 
is the integer of [4] —1 for the nth image. For convenience we have omitted the factor e™ 


which would otherwise appear in (2-11) and (2-12). 

Two special cases exist in the upper half space. For the case of the direct path of the 
source to the receiver (n=1), R is 1. In the case of a single reflection off the surface and no 
reflections with the bottom (n=2), R is again 1. The phase inversion of the ray is computed by 
the (-1)"7") term of (2-11) and (2-12). 
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FO 


Since the model is frequency independent, a scaling distance X, is used. X, is defined as 
the distance where the lowest normal mode attains cutoff when the bottom is fast (c2>c}). 


This cutoff distance is calculated by 
Xx -1] € 
a a a <co'(5] ae 
2sin 8, tan B C, } 
so, 
X. = 7k sin @ tanB oy 
For a slow bottom an analogous scaling, convenient for the computer program is 
n -1{ ©) 
XxX, = —— 6, = 4 2-15 
MX. 2 tan 6, tanB poe (2} poe 
so, 
X= Theta 8 tanB (2-16) 
c 2k tan 6 tan 
For the total pressure at the receiver, the upper and lower contributions are combined, 
thus 


P= P,+P, (2-17) 


il 





Ii. DOUBLET FORMATION 


A. GENERAL DESCRIPTION 
An acoustic doublet is a pair of acoustic images 180° out of phase. The first doublet 
consists of the source and its reflection, the 2nd upper image. This is called a neutral 


doublet (n=0). The reflection coefficient is 1 because its midpoint does not encounter a 





Figure 3-1 - Doublet Image Structure 
boundary. The 3rd and 4th upper images are the Ist upper doublet (n=1). The 5th and 6th 


upper images are the 2nd upper doublet (n=2). The Ist and 2nd lower images are the ist 
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lower doublet (n=1) and the 3rd and 4th lower images are the 2nd lower doublet (n=2). 
This procedure is continued until the opposite surface is reached. See Figure 3-1. 


The total number of acoustic doublet pairs is given by 


mn( 22") -2 
N=Int ; (3-1) 





At worst, only three images will be missed and they will be in the higher order images of 
the source so their contribution is negligible due to all the reflection coefficients 


encountered [Ref. 25]. 
B. DEVELOPMENT 
The sound propagation for a spherically spreading point source is 


p= A giter-b) (3-2) 
r 


where A is the pressure amplitude, r is the range from the source and k is the wave 


number o/c [Ref. 11]. Wher: two sources are combined to form a doublet, the above 


equation can be expanded to 
pa At pitar-b) _ AL ,ior-k,) (3-3) 
r r 
where the subscripts “+” and “-” refer to the upper and lower sources of the doublet and 
r,=rt+Ar A =A+AA (3-4) 
r_=r-A4r A_=A-AA 


The delta values are the incremental difference between each source compared to the 


theoretical point source. Figure 3-2 illustrates the doublet close-up. The distance d is the 
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separation of the images of the doublet pair and r is the distance of the doublet midpoint to 


the receiver. 


P(r, 8, t) 





» 


Figure 3-2 - Doublet Close-Up 


The pressure equation can now be rewritten as 








A+AA A-AA 
— A itar-tr) A jer.| A - jkr, 
aed r-Ar r+ ar fF 2) 
r r 


1. Balanced Doublet 
The Ist upper image has an reflection image of equal amplitude and opposite 
phase which is the second upper image. In the far field, the image pair is treated as a single 
doublet source. See Figure 3-3. 
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Receiver 





Figure 3-3 - Geometry of Neutral Doublet 


The pressure equation for a balanced source can be approximated by 
P= ca a = erio®)| (3-6) 
r 
' The geometry of the doublet is given in Figure 3-2. Ar is 
Ar= <sin(6) = r,sin(y)sin(6) (3-7) 


By substituting Equation 3-7 into Equation 3-6 and factoring one gets 





P= 2 em elie ( eltrisin(y sin(@) _ gir ene) (3-8) 
By using the relationship of 
~e* 


where x=kr;sin(y)sin(6), the result is 
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P= 2jAsincKr, sin(y)sin(@))e"-”” 
r 


2. Unbalanced Doublet 
For the unbalanced doublet, the amplitude relationship is expressed as 


A At A,-A. 


‘0 AA= 
2 


(3-11) 


The unbalanced doublet can be viewed as a combination of a balanced doublet on which is 


superimposed a like-charge doublet [Ref. 28]. See Figure 3-4 


nth Image Plane 


Figure 3-4 - Geometry of Unbalanced Doublet 
The pressure equation is 


P= [2 j As sinter sin(Y)sin(@)) + 2S costkr sin(y)sin( ane (3-12) 


which can be rewritten as 








P=2j2 
r 


sinc sin(y)sin(@)) — j = 


o 


cos(kr, sin( rin (or-#) (3.13) 


3. Distance from Doublet to Receiver 
The distance from the doublet formation to the receiver can be calculated using the 


law of cosines. For a three dimensional wedge, the following equation is used 


re =r, +1} + Yo — 27,7, cos(2nB + 5) (3-14) 
where r; and rz are the ranges from the apex to source and receiver respectively, yg is the 
cross-slope range, B is the wedge angle, 5 is the source angle, and n is the doublet 
formation number. The "+" and "-" signs are reference to the upper or lower images in 


the doublet formation. See Figure 3-5. 





Figure 3-5 - Doublet Three Dimensional Geometry 
Equation 3-13 can be rewritten by adding and subtracting 27;r2. This yields the 


following expression 


re sri —2nr, +r) + ye + 2r,r, — 2r,r, cos(2nB t 5) (3-15) 


By combining terms and factoring, the following is the result 





2 
7p =(r,-7,) [14+ +—25% (1 -cos(2np + 6)) (3-16) 
(r, -r,) (r, -r,) 
or 
ye 2r,r. 
ry, = (Fy — Fal [+27 + + 5 (1- cos(2nB t 5)) (3-17) 
(r, -r,) (r,-7, 
4. Pressure From Each Doublet 
A generalization of Equation 3-13 can be made by substituting A, and AA, for Ag 

and AA where 

A= AeA A, = Aen (3-18) 
and 6, for 8 where 

6, =2nBtd (3-19) 


so the pressure from each doublet is 


P,, =2j A inc sin(y)sin(9,,))— j “Ais cos kr 150 On » fr 


at 4 


(3-20) 
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IV. QUADRUPLET FORMATION 


A. QUADRUPLET GENERAL DESCRIPTION 
A quadruplet is formed by combining the upper and lower doublets. The visualization 


Ce.) 
\e 3] 
Primary Doubiet eel) 7 





Figure 4-1 - Quadruplet Illustration 
of the quadruplet is illustrated in Figure 4-1. The primary doublet (the source and its first 
upper image) are considered separately. The first quadruplet consists of the third and 
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fourth upper image (first upper doublet) and the first and second lower image (first lower 
doublet). The second quadruplet consists of the fifth and sixth upper images (second upper 
doublet) and the third and fourth lower images (second lower doublet). The process is 
continued until the opposite surface is reached. 

The number of quadruplets formed is 


N= in{ 2) (4-1) 


At worst, only 3 images will be missed and they are highest order images so their 


contribution is negligible due to all the reflection coefficients encountered. 


B. QUADRUPLET DEVELOPMENT 
1. Balanced Quadruplet 
The complementary doublets are 180° out of phase. The pressure of the 
quadruplet is given by combining the nth upper and lower doublets from Equation. 3-10 of 
the balanced doublet to form 


P,=2jA2] sin(kr, sin(y)sin(6,,))e 
r, (4-2) 


—sin(ky, sin(7)sin(@,_))e*~ JeX*-”) 
By substituting Equation 3-19 for 6,4 and using the approximation that sin(y)=y, the 
balanced quadruplet pressure equation is 


P, =2jd[ sin(kr,ysin(2nB + 5)" 
r, (4-3) 


—sin(kr,ysin(2nB — 5))e"*™"~ Je“) 
2. Unbalanced Quadruplet 
The amplitudes of the upper and lower doublets are not equal because of the 
distance differential between them as seen in Figure 4-2. This leads to the concept of the 
unbalanced quadruplet. 


nth Image Plane 


nth image Plane Bo AB 





Figure 4-2 - Illustration of Unbalanced Quadruplet Amplitude 
The Ar terms are derived from the comparison of the distance from the primary 
doublet to the midpoint of the upper and lower doublets in each quadruplet (Figure 4-2). 
This ensures that all phase angles are calculated from a common reference point. Using 
Equation 3-17, by setting yg=0, we can make the following approximation, 


Tat =\r, —~V, 1+ aa —conns 2) (4-4) 
1 2 
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rae =r, -nfio tte —cos(2nB t 5)) 
172 


Now Ar, can be approximated using the following relation 


Ar, 


2 =e VN Bc a el 
mah 


Midpoint of Neutral 
Doublet 


Figure 4-3 - Illustration of Ar in a Quadruplet 


The cosine term in the above expression can be expanded using a trigonometric 
identi 
cos(2nB + 5) = cos(2nB)cos(5) + sin(5)sin(2nB) (4-7) 
Small angle approximations can be used for 5 where cos(5)=1 and sin(5)=5 so 


cos(2nB + 5) = cos(2nB) + Ssin(2nB) (4-8) 





Equation 4-8 can be substituted into Equation 4-6 which can be substituted into 


Equation 4-3 to yield 


nh 


P= Diy getty {24 \0-conaab)y sin(kry sin(2nB + 5))e7*{ 2 Jalan 
r 


. (4-9) 


Ah 
Y 


—sin(kny sin(2n B- 6)e*( 7 = Jésin(ana) } 
For convenience, the following terms are defined 
¢ = kru(1-cos(2nB)) 
d = kr 6sin(2nB) 


__h 


li -a 
Equation 4-9 now becomes 
p= 2 ih, ell-re-# | sin(knysin(2np + 6))e™ 
—sin(kjy sin(2np - 5))e™ | 


The amplitudes of the upper and lower images of the doublets are not equal. By separating 


(4-11) 


the amplitudes and distributing it into the main equation in a process similarly done in 
Equation 3-12 and 3-13, the unbalanced quadruplet equation is the result 


2) jtwab) = 
P= Leia tr) i 
r 


[sins anp +6))- j-Aecos( kr ysin(2nf + a) jem (4-12) 


- Bin rsin(2nB a 5)) =J -cos( kr :ysin(2nBp - 6 ym 


where A, and B, are the amplitudes of the upper and lower nth doublets at the midpoint 
and AA, and AB, are the differences in amplitude of the nth upper and lower doublet 


images. 





if 
: 
i 
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V. REFLECTION COEFFICIENTS 


A. REFLECTION COEFFICIENTS 
The reflection coefficient is discussed in Chapter 2 with the image theory 
development. The reflection coefficient from [Ref. 29] is 
P2 sin(6,,) ~ <sin(6,) 
c. 


= 4 (5-1) 
Pe sin(Q +sin(@ 
Pa sin(0,)+ 2sin(6,) 


It can be rewritten as 


sin(6, ) 





~ sin(@,,,) 
ae se) a 


sin(0,. ) 





where 6,, is the grazing angle of the nth image going through the mth plane and 


b= c= (5-3) 
Pp C 


Sin(@,) can be expanded and the reflection coefficient can be rewritten as 
yl -c’ cos?(8, 
“a ae : - 
Run = — (5-4) 
1-c’ cos?(6,,) 
be+ - 
sin(8,,,) 

The term under the radical can be manipulated and written in terms of the critical angle, 6,, 
from Equation 2-13 
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2 
i-c?cos(8,,,) = 1 re) = fin On) =si0() (5-5) 


7 cos’( 


Equation 5-5 can be combined into Equation 5-4 and simplified. 





; 2 
p—.[y—| Sinl@) 
sin(8,,, ) 
ee eraaraGcy 5-6) 
pes sin(8, ) 
sin(6,,,) 
By using the substitution 
_ sin(6,.,) ; 
ee sin(@, ) oD 


the reflection coefficient takes the form 


bit 
(5-8) 


B. SLOW BOTTOM REFLECTION COEFFICIENTS 
A slow bottom is when the speed of sound in the bottom material is slower than in the 


water (c,>c2). This is an easier case to handle because there is no critical angle constraint, 
the reflection coefficient is real and diminishes quickly to zero. 
Equation 5-8 can be expanded to by using 


2 
of c-l (5-9) 


sin(@,) =1- a3 





resulting in 
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———— 


1-c? 
c’ sin(6,,) 


Dividing the top and bottom by 1/b and making an approximation with the square root 


resuits in 


Equation 5-13 can be approximated by 


2be_» 
R,, = re” 


Substituting 








2bc 


ee 5-16 
Vice (5-16) 


The final product is 


Ry, =e (5-17) 
Each ray intersects the bottom 2n-1] times. The product of the reflection coefficients 


will be 
R, = []49.] (5-18) 


Expanding the equation yields 


2n-1 


R = ERG Serre re ee (5-19) 


n 
m=1,3,... 


where 7 =€+6. Multiplying exponentials with the same base is equivalent to add their 


exponents. The terms in can be rearranged to 


2a-l 
I] A(¢,.] = e-NaNg a[1+345...(20-1)]8 (5-20) 


m=1,3,... 
By noticing that the sum of the expontentials at a index is equal to the index squared, the 
result is 

Rs ene (5-21) 
For small angles, the first term in the exponential is very small compared to the second so 


the cumulative reflection coefficient for a slow bottom is 
R, =e" (5-22) 


C. FAST BOTTOM REFLECTION COEFFICIENTS 
A fast bottom is when the speed of sound in the bottom is faster than the speed of 


sound of water (c2>c,). At angles below grazing, the acoustic energy will propagate with 
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100% reflectance (Rym=1) until 6, (see Equation 2-13) is reached then will decay. See 


Figure 5-1. 





v 





Fast Bottom Reflection Coefficient Curve for 3 Deg Wedge 


Abs(Reflection Coefficient) 





5 10 15 20 25 
Nth Quadruplet 





Figure 5-1 - Typical Fast Bottom Cumulative Reflection Coefficient Curve 


The development follows until Equation 5-8. By reversing the order of terms 


(5-23) 


(5-24) 





Using the Maclaurin series [Ref. 30] of 
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In a = 2 tanh"'(z) 


Equation 5-22 results in 


By taking out the 1/x? term and inverting the tanh"! 


ee Gel 


Simplifying the above equation results in 


ee 


-2j wo Ld , ) 
Run = ce (5-29) 
The cumulative reflection coefficient is the cumulative product of all previous reflection 


coefficients and will take the form of 


R, = i Rn Es al fs) 


a=] n=] 





VI. GRAZING AND RECEIVER ANGLES IN QUADRUPLET 


A. RECEIVER ANGLE 


The receiver angle is calculated using the law of sines. The relation can be seen in 


Figure 6-1 
nh r nt 
a aE tae (6-1) 
sin(é,.)  sin(2nB +d) 
By solving for sin(€,) and using the small angle approximation 
€,, =—Lsin(2nB + 6) (6-2) 


Tat 
Likewise the receiver angle for each image of the doublets can be solved for. The receiver 
angle for the upper half-space is given by 


at = eee; Se < (6-3) 
sin(é,,.)  sin(2nB+65+ 7) 
sin(¢,,.)=—1-sin(2np + 5+ ”) (6-4) 


n+t 


Using a the sine expansion of the sum of two angles, the above equation results in 
sin(e,,,) = <1 [sin(2np + 5)cos(y) + cos(2nB + 5)sin(7)] (6-5) 
Le 


Making the assumption cos(y)=1, sin(y)=y, sin(€)=€, r,,, =7,, and substituting sin(E,) 


from Equation 6-2 
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é,,,=&, + — ycos(2np + 6) (6-6) 


nt 


for the upper half space and 


E,- 


.=e,t a ycos(2nB - 5) (6-7) 


for the lower half space. 


B. ANGLE OF INCIDENCE FOR UPPER AND LOWER DOUBLET 
IMAGES 


The angle from the surface to each virtual bottom is 


6,, = (2n-1)B (6-8) 


The angle of incidence of the doublet midpoint is 


6,, =(2n-1)Btdt+e (6-9) 
The angles of incidence with the virtual bottom is different with each upper and lower 
doublet image. The angles of incidence for the upper half space are given by 


8,44 + (2n , )B +5+ Enat (6-10) 
By substituting Equation 6-5 into Equation 6-10 and making the standard small angle 


assumptions 
6... =(2n—-1B+5+ = [sin(2np + 5)(1)+ ycos(2nB + 5)] (6-11) 
Expanding the sine term and making the standard small angle assumptions again 
6,44 =(2n-1)B +5+ zi [sin(2nB)ay) + 5cos(2nB) + ycos(2nB + 5)|(6-12) 


Rearranging the terms result in 
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6,,, =(2n~1)B + <t sin(2nB) + 41 es - cos(an) 


(6-13) 
ot ycos(2nB + 6) 


a+ 


An assumption can be made that if r, << r,,, the sine term in the above equation is small 


so the it can be neglected. The result of Equation 6-13 is 


6,,.=9,,+ 41 + cos(anb) + — ycos(2nB + 5) 


n+ n 


The equivalent equation for the lower half space is 
6,5 = 8,, + i - “cotanb) $41 ycos(2nB — 6) 
r., r, 


Expanding the cosine term results in 


n+ 


6,44 =9,,+ a + “-cox(an)| 


(6-16) 
a4 y|cos(2nB)cos(5) — sin(2nB)sin(5)| 


Using the standard small angle approximations and assuming the sine term is small 
0,,, = 9,, + a1 + “i cos(an6) + = ycos(2nf) (6-17) 
for the upper half-space and 
0, =9,,+ a + 4 cos(anb + = ycos(2nB) 


for the lower half-space. 


C. AMPLITUDES 
The amplitude of each image in the quadruplet needs to be calculated due to their 
differences. The coefficients of Equation 6-17 and 6-18 are used where 





a, =1+-Lcos(2nB) 
r, 


b, = -cos(2nB) 
‘, 
The reflection coefficient can now be written for each image 


R,, =R,e"” (6-21) 


expanding the above equation and using the sine approximation for an angle 


R,, = Re e*"™™ (6-22) 


Each of the image amplitudes can now be expressed 
A, = R ennmd ener = R, Pend enmnr 


A. es R, end enmaY B., = R, end ear 
Using Equation 3-11 and the definition of 


e+e" 


: e*~e”* 
inh(x) = “$—*— 
2 IE 


cosh(x) = 
where x=nab,6 gives 

AA, =R,sinh(nab,d)e"™ AB, = R, sinh(nab,d)e"™" 

A, = R, cosh(nab,5)e"™-” B, = R, cosh(nab,5)e"™* 
The equations combine to form 


4 mtaah(n,7) pT tb(nab,7) 


Now the full quadruplet equation can be written as 





N * ; 
P=P+ >7R, cosh(nas, yee # 


7 haan 8 id (sin(kr,ysin(2nB + 5))—- jtanh(nab, y)cos(kr, ysin(2nB + 5))) 
—ea-5+ie (sin(kr, ysin(2np — 5)) — jtanh(nab, y)cos(kr, ysin(2nB — 6))) 


(6-27) 


where P, = “L sin (kr 76). 





VII. PROCEDURE AND RESULTS 


This thesis had three objectives: 
1. Convert the WEDGE/URTEXT to MATLAB™ and test it in all cases (slow and 
fast bottom, cross-slope, and with bottomloss). 


2. Linearize the quadruplet expansion code by Joyce (Ref. 27] 


3. Extend the quadruplet expansion to the fast bottom case with a approximation of 
the Rayleigh reflection coefficient. 


A. CONVERT WEDGE/URTEXT TO MATLAB 

The image model baseline program is called WEDGE. It is written in GWBASIC by 
Professor A. B. Coppens and is run on a IBM XT(8088) desktop computer. WEDGE 
was translated into FORTRAN by Nassopoulos [Ref. 25]. It was renamed URTEXT and 
is on the mainframe computer at the Naval Postgraduate School. Translating it to a 
language that is fast and can be run on a high performance desktop computer or 
workstation was necessary to make WEDGE data readily available for comparison to the 
quadruplet expansion. MATLAB™ [Ref. 31] was chosen because it is a high level 
scripting language which runs on a variety of computer operating systems including VAX, 
MS-DOS, UNIX, and Macintosh and that it is the mathematics program that the author is 
most familiar with. 

Validation of the MATLAB™ version of WEDGE, called WEDGEMAT, was 
-amducted by comparing the results of WEDGE and URTEXT to those of 
WEDGEMAT. 

Tables 7-1, through 7-4 show the results of the comparison of URTEXT with 
WEDGEMATt. WEDGEMATrr is the a version of WEDGE that ends the calculation 
process when the pressure amplitude is less than 0.00000001. The major error at the 
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surface is due to the difference in precision of the different languages. FORTRAN uses 


single precision variables, accurate to six decimal places. MATLAB™ is precise to 13 


places, therefore at the small pressures near the surface the round-off errors are extreme. It 
should be noted that the results of Nassoloulos are normalized pressures where the actual 
pressure is multiplied by the source scaling distance rj. 

Tables 7-5 through 7-8 show the comparison of WEDGE to WEDGEMATtr. The 
WEDGE program in GWBASIC uses variables precise to 12 decimal places so it 
compares more favorably with the MATLAB version. Runs were done with fast bottoms, 
bottomloss, and cross-slope to test accuracy. The results are very favorable. 

A run was made on an 10° slope with the data taken from the surface to bottom 
sampled every 0.05°. Figure 7-1 shows pressure falls off linearly as the sound nears the 
surface (bottom of graph). This implies the model is valid at the surface. 


N 


Receiver angle measured 
from the surface 


Receiver Angle (Deg) 
in 


—= 


Figure 7-1 - Pressure vs. Receiver Angle 
WEDGEMATrr has undergone basic proofing. When compared with previous 
versions in FORTRAN and GWBASIC, it is more flexible both in platforms run, model 








variables, and speed. Most importantly, now it can be used to check the quadruplet 


expansion data to verify results. 


TABLE 7-1 SLOW BOTTOM COMPARISON OF URTEXT AND 
WEDGEMATTR 
B=3  y=1 ry=l m=10 yo=0 pi/p2=0.9 
c1=1500m/s c2=1475m/s a/k2=0 


Pressure Amplitude Phase Angle 


Absolute 


Difference 


; 0.002471 . 000004 . , : 
pd press ite noone | Agree Dit is _| 


Note: Normalized pressure equation used where P = 7,./P? + P? 
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TABLE 7-2 SLOW BOTTOM COMPARISON OF URTEXT AND 
WEDGEMATTR 
B=4  y=1l  orys5  12=12 yo=0 p1/p2=0.7 
c1=1500m/s c2=1420m/s oa/k2=0 : 


: ; opens 
Pressure Amplitude Phase Angle 


W Receiver Absolute Absolute 


Angle Difference MATLAB _ Difference | 


0.040927 0.041003 
0.394236 0.394324 
0.701851 0.701860 
0.873027 0.873049 
0.915180 0.915294 
0.887048 0.887038 
0.851987 0.851901 
0.810276 0.810180 


0.704744 0.704662 





Note: Normalized pressure equation used where P = r,./P? + P? 
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TABLE 7-3 SLOW BOTTOM COMPARISON OF URTEXT AND 
WEDGEMATTR 
p=5 y=4 r1=2 12=400 yo=0 p1/p2=0.6 
c1=1500m/s c2=1410m/s o/k2=0 


Pressure Amplitude Phase Angle 


| 
| 
| 

























A Receiver Absolute Absolute 












\ 
h Angle Difference | 







0.000087 


0.0008 18 









0.001532 







0.002151 







0.002537 







0.002704 


0.002639 






0.002276 


i . .000010 .. : . 
a, Avg Pressure Diff= 0.000009 Avg Phase Diff= 3.7" 


Note: Normalized pressure equation used where P = 1,./P? + P’ 






TABLE 7-4 SLOW BOTTOM COMPARISON OF URTEXT AND 
WEDGEMATTR 
B=6 y=2 1r1=0.8 r2=5 yo=0 p1/p2=0.8 
c1=1500m/s c2=1450m/s o/k2=0 


sii hadi 


Receiver Absolute Absolute 
Angle Difference Difference 
0.000372 0.000377 
0.003728 0.003726 


0.007214 0.007205 


| 
| 
| 
| 
| 


0.010212 0.010210 
0.012559 0.012557 
0.014134 0.014136 
0.014946 0.014947 
0.015168 0.015166 

0.015224 


0.015832 





Note: Normalized pressure equation used where P = 7,./ P? + P? 
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TABLE 7-5 FAST BOTTOM COMPARISON OF WEDGE AND 
WEDGEMATTR WITH BOTTOMLOSS 
p=3 y=l rm=l 12=10 yo=0 p1/p2=0.9 
c1=1500m/s c2=1666m/s o/k2=0.0001 ss 


Pressure Amplitude Phase Angle 


Absolute Absolute 





Difference Difference 


0.01255 


0.12377 


0.23681 


0.32982 


0.39272 


0.40585 


0.34832 


| 
0.42214 | 
| 
| 
| 


2. 0.14819__—O. 00003 | 
—. Avg Pressure Diff= 0. Avg Diff= 0.5" | 
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TABLE 7-6 FAST BOTTOM COMPARISON OF WEDGE AND 
WEDGEMATTR WITH BOTTOMLOSS 
p=3  y=1 r=l m=100 yo=0 21/p2=0.9 


c1=1500m/s c2=1875m/s_ o/k2=0.0001 


Pressure Amplitude Phase Angle , 


| 
Absolute Absolute | 


Difference Difference 


0.00199 


0.02062 


0.03966 


0.05535 
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TABLE 7-7 FAST BOTTOM COMPARISON OF WEDGE AND 
WEDGEMATTR IN CROSS-SLOPE WITH BOTTOMLOSS 
B=3 y=1 ry=l m=10 yo=4 p1/p2=0.9 
c1=1500m/s c2=1666m/s o/k2=0.0001 


Pressure Amplitude Phase Angle | 
| 


Absolute Absolute 
Difference Difference 
0.01024 0.01094 
0.10162 0.10227 
0.19699 0.19716 
0.27908 0.27931 


0.34418 0.34418 


| 0.39920 0.39896 


0.39954 
0.35799 
0.27659 


0.17036 


| 
[sata oz | 
Avg Pressure Diff= _ 0.00029 Avg Phase Diff=__0.2° | 





TABLE 7-8 FAST BOTTOM COMPARISON OF WEDGE AND 
WEDGEMATTR IN CROSS-SLOPE WITH BOTTOMLOSS 
B=3 y=0.5 ry=l1 m2=10 yo=8 pj/p2=0.9 
cj=1500m/s c2=1666m/s a/k2=0.0001 


Absolute 
Difference | 
0.00501 
0.04911 
0.09495 
0.13463 
0.16949 


0.18845 


0.18554 


0.16382 





B. LINEARIZATION OF QUADRUPLET EXPANSION PROGRAM 

The original quadruplet expansion program was written by Joyce [Ref. 27]. His 
program is used as a basis for the fast bottom case examined later in this chapter. The 
linearized program is in Appendix A-1 and A-2 and has incorporated into it the fast 


bottom, reflection coefficients of Equations 7-1 and 5-28. 
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The original quadruplet expansion code used "For Loops" to compute its products. 
The code was linearized to speed processing time. A speed factor of 8 was achieved with 
the same output results [Ref. 27}. Table 7-9 shows relative execution time for a 3° wedge. 
This table is only meant to give the reader a ‘feel’ of the time involved in running the 
programs. The WEDGE program in GWBASIC was run on a IBM XT (8088) with 
FPU. All other programs were run on a Mac IIsi (68030 @ 20 MHz) with Applied 


Engineering 32k RAM cache card (with 20 MHz FPU) using MATLAB™ 3,5, 
TABLE 7-9 EXECUTION TIME OF IMAGE PROGRAMS 





Another method to speed calculation is to use the "For Loop” and put in a "Break" 
command to stop the calculation process after the pressure for a certain image or image 
doublets, in the quadruplet case, falls below a specified level. This was done in the 
WEDGE (truncated version) and WEDGEMATItr. This especially useful in a slow 
bottom case where the first 3 of 30 quadruplets (for a 3° wedge) would be 98 percent of the 
amplitude. This idea was not explored for the quadruplet expansion case. 


C. FAST BOTTOM REFLECTION COEFFICIENT APPROXIMATION 
The Rayleigh reflection coefficient without absorption is 











pic sin(8,..) 
Ry =—— , (7-1) 
(=) cos’(6,,,)—1 
Por _ , VAM 
pc, sin(6,,) 


can be approximated by Equation 5-28 which is 


R= _ue| ad 
For a fast bottom 3° wedge with the following parameters 

c, =1500m/s p, =1Kg/ m? 

Cc, = 1666m/s p, =LAL11Kg/m? 
the reflection coefficients are displayed in Table 7-10. In a 3° wedge, 30 reflection 
coefficients would be calculated. A fast bottom requires the first 10 terms since the 
cumulative product is multiplied with the amplitude of the quadruplet. The cumulative 
product after the tenth term is very small and contributes little to the final pressure as can be 
seen in Figure 5-1. 
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TABLE 7-10 COMPARISON OF RAYLEIGH AND ARCTAN 
APPROXIMATION REFLECTION COEFFICIENTS 
Rayleigh Reflection Tan"! 
Coefficient Approximation 


— 


2 
3 
4 
5 
6 
7 
8 
9 


— 
So 


The correlation of the two reflection coefficients are good until the third quadruplet. 
Then the Tan-! approximation fall off much more quickly than the Rayleigh values. This 
can be seen in the complex plane plot of the reflection coefficients in Figure 7-2. The first 
five reflection coefficients are numbered in each case. 





Plot of Complex Reflection Coefficient 


x - Rayleigh 
* ~ ArcTan Approx. 


Figure 7-2 - Comparison of Complex Reflection Coefficients 
TABLE 7-11 COMPARISON OF THE CUMULATIVE RAYLEIGH AND 
ACTTAN APPROXIMATION REFLECTION COEFFICIENTS 


Rayleigh Reflection Coefficient 








The most important factor is the cumulative :¢ ‘ection coefficient which is shown in 






Table 7-11. There is a wide disparity between the actual and approximated reflection 





coefficients. The correlation is good until the second quadruplet, then the products diverge. 







: The reason for this is the inaccuracy of each grazing with the bottom is amplified by the 






complex multiplication of the cumulative product. 






A test was run comparing the amplitude and phase at the receiver. The output is 






displayed in Table 7-12. Data were taken from the surface to the bottom in 0.5° 


increments. The differences are attributed to the inaccuracy of the reflection coefficient 






after the second quadruplet. With the fast bottom, the equations are very sensitive to the 







phase of the reflection coefficients. Even though the reflection coefficients have an 





amplitude of 1, the phase can cause constructive and destructive interference. Also with a 






fast bottom, many of the terms in the quadruplet expansion are complex there by 
magnifying the effect of the phase inaccuracy. 
TABLE 7-12 COMPARISON OF AMPLITUDE AND PHASE OF USING 
RAYLEIGH AND ARCTAN APPROXIMATION 















Receiver Angle 
(from Surface) 


Rayleigh Reflection Coefficient 








fpr rR 


Now one still has to compare the results of WEDGEMATrtr, to the quadruplet 
expansion. The quadruplet expansion with the Rayleigh reflection coefficient was tested in 
the same conditions as in the above table. The results are displayed in Table 7-13. 

TABLE 7-13 COMPARISON OF AMPLITUDE AND PHASE OF 
QUADRUPLET EXPANSION USING RAYLEIGH AND IMAGE MODEL 


Receiver Angle | Rayleigh Reflection Coefficient WEDGEMATIr 





The comparison of the image model (WEDGEMATKtr) and the quadruplet expansion 
using the Rayleigh reflection coefficient reveals an error factor of about 3 in amplitude. 
This can be attributed to the many approximations made in deriving the equation for the 
quadrupiet expansion which was covered in Chapter III to VI of this thesis. Many of the 
approximations are first order and are not good enough to get accurate results. 

The quadrupiet expansion using the exponential approximation for a slow bottom 
(Ref. 27] is limited to small wedge angle, the source close to the surface, and the receiver 
must be in the ‘far field’. The reason for the limitation to a small wedge angle (B~3) is to 
minimize the errors in the small angle approximation ( cos(y)=1 and sin(5)=5). The 
source needs to be close to the surface (upper half of the wedge) so the distance between 
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the upper and lower doublet pairs are reasonably close. The ‘far fie'd' requires r, >> 1r,, 


which is about 50 dump distances (X,). The reason for this is to insure the receiver angle 
(€) is small and the difference between the upper and lower images in the doublet are very 
small. These requirements are needed to insure that in the slow bottom case, where the 
reflection coefficient is real, and the cumulative product of the reflection coefficient is ~O by 
the third quadruplet (of 30) in a 3° wedge. In the fast bottom case, when the terms are 
complex and phase dependent, the cumulative reflection coefficient will go to ~0 by the 
ninth or tenth quadruplet. in this latter case there are many more complex operations being 
executed thereby magnifying the need for accurate approximations of the reflection 
coefficients and the terms in the quadruplet expression (Equation 6-27). 





Vill. CONCLUSION AND RECOMMENDATIONS 


Image theory is an ideal method for calculating the transmission loss in a shallow 
water (wedge shaped ocean) environment. It can be used in cross-slope, at all frequencies 
and in transitional cut off regions that are out of bounds for normal mode theories. 

This thesis has successfully converted the GWBASIC and FORTRAN code of 
WEDGE and URTEXT to MATLAB™, which is a high level scripting language. The 
original WEDGE program is accurate in fast and slow bottoms, cross-slope, and with 
bottomloss. The WEDGEMAT program retains all the same qualities but is much faster 
to run. The quadruplet expansion program for a slow bottom by Joyce (Ref. 27] was 
linearized to increase processing speed by a factor of 8 while retaining the same results. 

This thesis also incorporated a Tan~! approximation of the reflection coefficient for a 
fast bottom into the quadruplet expansion. Due to the inaccuracy of the reflection 
coefficient after the second quadruplet, the results were not favorable. It was also 
discovered that even with an accurate Rayleigh reflection coefficient, the first order 
approximations made in developing the quadruplet expansion equation (Equation 6-27) are 
not accurate enough for the fast bottom case. 

The next steps in this area of research is to find a better approximation for the fast 
bottom reflection coefficient and to increase the accuracy of the terms of the quadruplet 


expansion equation (Equation 6-27). 
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APPENDIX A-1 


QUADATAN 
Linearized Quadruplet Expansion of Method 
of Images. Based on PQUAD.m by Joyce. 


This program will provide pressure amplitude and 

the phase angle at the reciever. All distances are 
ratios of dump distances Xc. All angles are inputed 
in degrees. This model includes Rayleigh Reflection 
coefficient for fast bottom case and exponential 
approximation for slow bottom case. 


By Pat Takamiya 


dP dP OP GP dP dP dP dP dP dP dP dP UP dP 


clear; cle; clg; 

% Input Parameters 

B=input ('Enter bottom wedge angle in deg (beta) ='); 

G=input (‘Enter source angle from surface in deg (gamma) ='); 
D=input('Enter receiver angle from surface in deg (delta)='); 
di=input('Enter Density ratio, water to bottom ='); 

CC=input ('Enter sound speed ratio, water to bottom ='); 
rl=input('Enter range of source from apex (rl) ='); 
r2=input('Enter range of receiver from apex (r2) ='); 

yO=0; 

% 

% Determine Number of Image Quadruplets 

% Convert Angles from Deg to Rads 

N1=£1x(90/B) ; 

B=B* (pi/180) ; 


=G* (pi/180); 
D=D* (pi/180); 
% 
% Determine Scaling for Fast or Slow Bottoms 
tqq=tan (B) ; 
if cc<1 % Fast Bottom 


tqql=acos (CC) ; 
tqq2=sin(tqql) ; 


else 
tqql=acos (1/CC); % Slow Bottom 
tqq2=tan (tqq1) ; 
end 
tqq3=2*tqq2*tqq; 
t4=pi/tqq3; % Scaled wave number 
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% 

% Calculate Constants and Primary Doublet 
% 

alpha=2* (1/d1) / (sqrt ((CC*2)-1)); 

r=sqrt ((r1°2) +(r2*2) + (y0%2) - (2*r1l*r2*cos(D))); 
r3=abs(r2-r1); 

r4=r1*r2/(r3%2); 

mu=r2/r3; 

q=t4*r1*G*D; 

pl=(2/r)*sin(q) ; 

£=0; 

% 

% Quadruplet Summation 

% 


n=(1:N1); 

s(n)=n; 
th(n)=2*n*B; 
a(n) =t4*r1*D*sin(th(n)); 
phi (n) =t4*r1*mu* (1-cos(th(n))); 

% Reflection Coeff 
ru(n) =sqrt ((r1*2) + (r2*2) + (y0%2) -2*r1*r2*cos(th({(n)+D)); 
r1(n)=sqrt ( (r1*2)+(r2*2)+(y0%*%2) -2*r1*r2*cos(th(n)-D)); 
r5(n)=r1./ru(nj; 
r6(n)=r1./r1(n); 
al(n)=1¢r5(n).*cos(th(n)); 
bl1(n)=r6(n).*cos(th(n)); 
if ccC>l1 
rr (n)=(B*alpha.*(n.%*2)); 
R(n)=-exp(-rr(n)); 
else 
re=(1/CC) *1/dl1; 
ec1=1/CC; 
gi (n)=sqrt (ccl*ccel*cos(th(n))-1); 
g2(n)=sin(th(n)); 
g3 (n)=j*g1(n) ./g2(n); 
R1 (n) =(re+g3 (n)) ./ (re~-g3 (n)); 
R(n)=cumprod({-R1(n)); 
end 


a(n) =n.*alpha*D.*al(n); 
b(n) =n. *alpha*G.*bl1(n); 


% Pressure 
ul(n)=cosh(b(n)); 
u2 (n) =exp(-j*phi(n)); 
u3 (n)=(-2/r).*R(n) .*ul(n) .*u2(n); 





% Upper Image 
vi (n)=sin(th(n)+D); 
v2 (n) =t4*r1*G*v1(n) ; 
v3 (n) =sin(v2(n)); 
v4(n)=cos (v2 (n)); 
v5 (n) =-a(n) -j*mu*d(n) ; 
v6 (n) =exp(v5(n) ); 
v7 (n)=tanh(b(n)); 
v8 (n) =v6(n) .* (v3 (n) -j*v7(n) .*v4(n)); 


% Lower Image 
wl (n)=sin(th(n)-D); 
w2 (n) =t4*r1*G*wil (n) ; 
w3 (n) =sin(w2(n)); 
w4 (n) =cos (w2 (n)); 
w5 (n) =a(n)+j*mu*d(n) ; 
w6 (n) =exp (w5 (n) ); 
w7 (n) =w6(n).* (w3(n)-j.*v7(n).*w4(n)); 
% 
% Summation 
p(n) =u3 (n) .* (v8 (n)-w7(n)); 
p2=sum(p(n)); 
% 


% 

ps=p1+p2; 

phiz=atan2 (imag(ps) , real (ps) ) ; 
phase=180*ph1/pi 

P=sqrt (real(ps)*2 + imag(ps) *2) 
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APPENDIX A-2 


% 

% QUADATAN 

% Linearized Quadruplet Expansion of Method 

% of Images. Based on PQUAD.m by Joyce. 

% 

% This program will provide pressure amplitude and 

% the phase angle at the reciever. All distances are 
% ratios of dump distances Xc. All angles are inputed 
% in degrees. This model includes ArcTan Approximation 
% Reflection coefficient for fast bottom case and 
exponential 

% approximation for slow bottom case. 

% 

% By Pat Takamiya 

% 


clear; cle; clg; 

% Input Parameters 

B=input (‘Enter bottom wedge angle in deg (beta) ='); 
G=input('Enter source angle from surface in deg (gamma) ='); 
D=input ('Enter receiver angle from surface in deg (delta) 
='); 

di=input('Enter Density ratio, water to bottom ='); 
CC=input('Enter sound speed ratio, water to bottom ='); 


rl=input('Enter range of source from apex (rl) ='); 
r2=input('Enter range of receiver from apex (r2) ='); 
y0=0; 

% 


% Determine Number of Image Quadruplets 
% Convert Angles from Deg to Rads 
N1=fix(90/B) ; 

B=B* (pi/180) ; 

G=G* (pi/180) ; 

D=D* (pi/180) ; 


% 
% Determine Scaling for Fast or Slow Bottoms 
tqq=tan (B) ; 
if cce<1 % Fast Bottom 
tqql=acos (CC); 
tqq2=sin (tqql) ; 
else 


tqqli=acos (1/CC) ; % Slow Bottom 
tqq2=tan (taq1) ; 
end 
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tqq3=2*taqq2*taq; 
t4=pi/tqq3; % Scaled wave number 
% 


% Calculate Constants and Primary Doublet 
% 
alpha=2* (1/d1) / (sqrt ((CC*2)-1)); 
r=sqrt ( (r1*2)+(r2°2) +(y0*%2) -(2*r1*r2*cos(D))); 
r3=abs (r2-r1); 
r4=rl1*r2/(r3%2); 
mu=r2/r3; 
q=t4*r1*G*D; 
pl=(2/r)*sin(q); 
£=0; 
% 
% Quadruplet Summation 
% 
n=(1:N1); 
s(n)=n; 
th (n) =2*n*B; 
a(n) =t4*r1*D*sin(th(n)); 
phi (n) =t4*r1*mu* (1-cos(th(n))); 
% Reflection Coeff 
ru(n)=sqrt ((r1%*2) + (r2%2) + (y0*2) -2*r1l*r2*cos(th(n)+D)); 
r1(n)=sqrt ( (r1*2)+(r2%2) + (y0%2) -2*r1*r2*cos (th(n)-D)); 
r5(n)=r1./ru(n); 
r6(n)=r1./rl1(n); 
al (n)=1+r5(n).*cos(th(n)); 
b1(n)=r6(n).*cos(th(n)); 
if cc>l 
rr (n)=(B*alpha.*(n.%*2)); 
R(n)=-exp(~-rr(n)); 
else 
x(n)=sin(th(n))/sin(acos(CC)); 
R1 (n) =-exp (-2*j*atan(((1/d1) .*x(n)) ./sqrt(1-x(n) .*2))); 
R(n) =cumprod (Ri (n) ) ; 
end 
a(n)=n.*alpha*D.*al(n); 
b(n) =n. *alpha*G. *bl1 (n) ; 


% 
% Pressure 
ui(n)=cosh(b(n)); 
u2 (n) =exp(-j*phi(n)); 
u3 (n)=(-2/r) .*R(n) .*ul(n) .*u2(n); 
% 


% Upper Image 
v1(n)=sin(th(n)+D) ; 
v2 (n) =t4*r1*G*vl1(n); 
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v3 (n)=sin(v2(n)); 

v4(n)=cos (v2 (n)); 

v5 (n) =-a(n) -j *mu*d(n) ; 

v6 (n) =exp (v5 (n)); 

v7 (n) =tanh(b(n)); 

+8 (n)=v6(n) .* (v3 (n) -j*v7(n) .*v4(n)); 


% Lower Image 

wi (n)=sin(th(n)-D); 

w2 (n) =t4*r1*G*wl (n) ; 

w3 (n) =sin(w2(n)); 

w4 (n) =cos (w2 (n)); 

w5 (n) =a(n)+j*mu*d(n) ; 

w6 (n) =exp (w5 (n) ); 

w7 (n) =w6 (n) . * (w3(n)-j.*v7 (n) .*w4(n)); 
% 
% Summation 
p(n) =u3 (n) . * (v8 (n) -w7 (n)); 
p2=sum(p(n)); 
% 
% 
ps=pl+p2; 
phi=atan2 (imag (ps), real(ps)); 
phase=180*ph1/pi 
P=sqrt (real(ps)*2 + imag(ps) “2) 
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APPENDIX B-1 


WEDGE for MATLAB 


Adapted from the BASIC version of WEDGE 
by Prof A.B. Coppens and URTEXT by G. Nassopolous 
pz Not Normalized 


This program will provide pressure amplitude and 

the phase angle at the reciever. All distances are 
ratios of dump distances Xc. All angles are inputed 
in degrees. This model includes cross~-slope and 
bottomloss terms. 


By Pat Takamiya 


% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 


clear; clc; clg; 
x 2 
B=input('Enter bottom wedge angle (beta)= '); 
G=input('Enter source angle from the surface (gamma)= '); 
D=input('Enter receiver angle from the surface (delta)= '); 
di=input('Enter Density ratio, water to bottom= '); 


CC=input('Enter speed ratio, water to bottom= '); 
ri=input('Enter range of source from apex (r1)= '); 
r2=input('Enter range of receiver from apex (r2)= '); 
BL=input ('Enter bottom loss coefficient (alpha/k2)= '); 
yO=input ('Enter cross-slope range (Yo)= '); 
% 
% Determine Number of Image Pairs 
% Convert Angles from Deg to Rad 
% 
N1=f£ix(180/B) ; 
B=B*pi/180; 
G=B-G*pi/180; 
D=B-D*pi/180; 
% 
% Determine Scaling for Fast or Slow Bottom 
% 
tqq=tan(B); 
if cc<l1, % Fast Bottom 
tqql=acos (CC); 
tqq2=sin (tqql) ; 
else 





tqql=acos(1/CC); % Slow Bottom 
tqq2=tan (tqql) ; 
end; 
% 
% Determine Constants 
% 
tqq3=2*tqq2*tqq; 
t4=pi/tqq3; 
c2=CC*%2; 
d2=(r1*2)+(r2*2)+(y0%2); 
r3=2*r1*r2; 
ql=1/sqrt (2); 
% 
% Define Lengths of arrays 
i=1:fix(N1/2); 
c(i)=zeros (1:£1x(N1/2) ); 
e(i)=zeros (1:£ix(N1/2) ); 
£(i)=zeros (1: fix(N1/2)); 
% 
% Determine Range to Receiver for each Image 
% 
n=1:Nl; 
n1=2:2:N1; 
n2=1:2:N1; 
ti(n1)=(n1*B) -G; 
t1(n2)=((n2-1) *B) +G; 
r8(n)=sqrt (d2-r3*cos(t1l(n)-D)); % 
r9 (n)=sqrt (d2-r3*cos(t1l(n)+D)); % 


3-D Law of Cosines 
3-D Law of Cosines 


% 
pl=0;p2=0; 


for n=1:1:N1, 
s2=(-1)*(fix(n/2)); 


% 
% Calculate Reflection Coefficient for Upper Space 
% 
wl=2*c2*BL; 
i1=f1ix((n-1)/2); 
for i=1:1:i1, 
s(i)=(r1*sin(t1(n)-2*i*B) +r2*sin(2*i*B-D) )/r8(n); 
if s(i)>=1, 
s(i)=1; 
end; 
c(i)=sqrt (1.0-s(i)%*2); 
t=s(i)/d1; 
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w0=-c2+c (i) %*2; 
y=sqrt ( (w0*2)+(wl%2)); 
z=(w0); 
if y<=z, 
y=2; 
end; 
yl=qi*sqrt (y+w0) ; 
y2=-ql1*sqrt (y-w0) ; 
zl=t-y2; 
z2=-yl; 
Z3=21/ ((z1*2)+(z2%2)); 
z4=-22/ ((z1*2)+(z2%2)); 
zl=tty2; 
z2=yl; 
z5=(z1*z3) -(z2*z4) ; 
z6=(21*24) +(z2*z3); 
e(i)=z5; 
£(i)=z6; 
end; 
% 
% Determine Pressure Contribution for the Upper Images 
% 
21=0; z2=0; 23=0; 24=0; 25=1; 26=0; 
if n>2, 
for i=1:1:i1, 
zi=e(i); 
z2=f (i); 
Z3=25; 
z4=z6; 
z5=(21*Z3) -(z2*z4); 
z26=(zZ1*z4)+(z2*z3); 
end; 
end; 
z1=z5; 
z2=26; 
t=t4*r8(n); 
z3=cos (t); 
z4=~-sin(t); 
25=(z1*z3) -(22*2z4) ; 
z6=(z1*z4) +(2z2*z3); 
pl=pl+ (s2*z5/r8(n)); 
p2=p2+ (s2*z6/r8(n)); 
il=i1+1; 
% 
% Determine Reflection Coefficient for the Lower Images 
% 


for i=1:1:i1, 








s(i)=(ri*sin(tl1(n) -2* (i1) *B)+r2*sin(2* (i1) *B+D) )/r9(n); 
if s(i)>1, 
s(i)=1; 
end; 
e(i)=sqrt(1.0-s(i)%*2); 
t=s(i)/dl; 
w0=-c2+c (i) %2; 
y=sqrt ( (w0*2)+(w1%*2)); 
z=(w0); 
if y<=z, 
Y=2; 
end; 
yl=ql* sqrt (y+w0) ; 
y2=-ql*sqrt (y-w0) ; 
z1l=t-y2; 
z2=-yl; 
23=Z1/ ((z1*2)+(22%2)); 
24=-22/ ((z1*°2)+(2z2%*2)); 
z1l=t+y2; 
z2=yl; 
z5=(z1*23)-(z2*z4); 
26=(z1*z4)+(z2*z3); 
e(i)=z5; 
£(i)=z6; 
end; 
% 
% Find Pressure Contribution for Lower Images 
% 
21=0;22=0; 23=0; 24=0; z5=1; z6=0; 
for i=1:1:i1, 
zl=e(i); 
z2=f (i); 
z3=z5; 
z4=26; 
Z5=(Z1*zZ3) -(z22*24) ; 
z6=(zZ2*Z3)+(z1*24) ; 
end; 
zl=z5; 
z2=Z6; 
t=t4*r9(n); 
z3=cos(t); 
z4=-sin(t); 
z5=(z1*z3) -(z2*z4) ; 
z6=(22*23)+(z1*z4) ; 
pl=p1+ (s2*z5/r9(n)); 
p2=p2+ (s2*z6/r9(n) ); 
end; 
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% 
% Total Pressure at the Reciever 


% 

tS=sqrt (p1*2+p2%2) % Pressure at Reciever 
% Phase Angle at the Reciever 

phaseang=atan2 (p2,p1)*180/pi % Phase Angle 


APPENDIX B-2 


WEDGE for MATLAB 


Adapted from the BASIC version of WEDGE 
by Prof A.B. Coppens and URTEXT by G. Nassopolous 
pz Not Normalized 


This program will provide pressure amplitude and 

the phase angle at the reciever. All distances are 
ratios of dump distances Xc. All angles are inputted 
in degrees. This model includes cross-slope, 
bottomloss, and is truncated at pz<0.00000001. 


By Pat Takamiya 


dP dP dP dP dP dP dP dP dP dP dP dP dP dP dP 


clear; clc; clg; 
% 
B=input('Enter bottom wedge angle (beta)= '); 
G=input ('Enter source angle from the surface (gamma)= '); 
D=input ('Enter receiver angle from the surface (delta)= '); 
dli=input('Enter Density ratio, water to bottom= '); 
CC=input('Enter speed ratio, water to bottom= '); 
rli=input('Enter range of source from apex (r1)= '); 
r2=input ('Enter range of receiver from apex (r2)= '); 
BL=input('Enter bottom loss coefficient (alpha/k2)= '); 
y0O=input (‘Enter cross-slope range (Yo)= '); 
% 
% Determine Number of Image Pairs 
% Convert Angles from Deg to Rad 
% 
N1i=f£1ix(180/B) ; 
B=B*pi/180; 
G=B-G*pi/180; 
D=B-D*pi/180; 
% 
% Determine Scaling for Fast or Slow Bottom 
% 
tqq=tan (B) ; 
if Ce<l1, % Fast Bottom 
tqql=acos (CC); 
tqq2=sin(tqq1) ; 
else 
tqql=acos (1/CC) ; % Slow Bottom 
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tqq2=tan (tqql); 
end; 

% 

% Determine Constants 

% 

tqq3=2*tqq2*taqq; 

t4=pi/tqq3; 

c2=CC*2; 

d2=(ri*2)+(r2%*2)+(y0%2) ; 

r3=2*xr1*r2; 

qi=1/sqrt (2); 

gqa=1; % Truncation Tester 
% 

% Define Lengths of arrays 
1=1:fix(N1/2); 
c(i)=zeros(1:£ix(N1/2)); 
e(i)=zeros(1:fix(N1/2)); 
£(i)=zeros(1:fix(N1/2)); 






















% 
% Determine Range to Receiver for each Image 


% 





n=1:N1; 
n1=2:2:N1; 
n2=1:2:N1; 
ti(n1)=(n1*B) -G; 
t1(n2)=((n2-1) *B) +G; 
r8(n)=sqrt (d2-r3*cos(t1(n)-D)); % 3-D Law of Cosines 
r9(n)=sqrt (d2-r3*cos(t1(n)+D)); % 3-D Law of Cosines 











% 
p1=0;p2=0; 










for n=1:1:N1, 
s2=(-1)*(f£ix(n/2)); 










% 
% Calculate Reflection Coefficient for Upper Space 
% 
wl=2*c2*BL; 
i1=fix((n-1)/2); 
for i=1:1:il, 
s(i)=(r1*sin (tl (n) -2*i*B) +r2*sin(2*i*B-D) )/r8(n); 
if s(i)>=1, 
s(i)=1; 
end; 
e(i)=sqrt (1.0-s(i)%*2); 
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t=s(i)/dl; 
w0=-c2+c (i) %*2; 
=sqrt ( (w0*2)+(w1%2) ); 

z=(w0); 

if y<=z, 

y=2Z; 

end; 
yl=ql*sqrt (y+w0) ; 
y2=-ql*sqrt (y-w0) ; 
z1l=t-y2; 
z2=-yl; 
Z3=Z1/ ((z1*%2)+(z2%*2)); 
Z4=-22/ ((z1°2)+(z2%2)):; 
zl=t+y2; 
z2=yl; 
z5=(z1*z3) -(z2*z4); 
Z26=(zZ1*z4)+(z2*z3); 
e(i)=z5; 
£(i)=z6; 

end; 


% 
% Determine Pressure Contribution for the Upper Images 
% 

Z21=0; 22=0; z3=0; 24=0; 25=1; z6=0; 


if n>2, 
for i=1:1:i1, 
zl=e(i); 
z2=£ (i); 
z23=2z5; 
z4=2z6; 
z5=(z1*z3)-(22*z4) ; 
z6=(z1*z24)+(z2*z3); 
end; 
end; 
z1l=z5; 
z2=2z6; 


t=t4*r8(n); 
z3=cos(t); 
z4=-sin(t); 
25= (Z1*z3) -(z2*24) ; 
z26= (Z1*z4) +(z2*z3) ; 
pl=pl1+ (s2*z5/r8(n)); 
p2=p2+ (s2*z6/r8(n)); 
i1=i1+1; 
% 
% Determine Reflection Coefficient for the Lower Images 
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for i=1:1:i1, 
s(i)=(r1l*sin(t1i(n)2* (il) *B)+r2*sin(2* (il) *B+D))/r9(n); 
if s(i)>1, 
s(i)=1; 
end; 
c(i)=sqrt(1.0-s(i)%*2); 
t=s(i)/dl; 
w0=-c2+c(i)%*2; 
y=sqrt ( (w0*2) +(w1%2)); 
z= (w0); 
if y<=z, 
y=2; 
end; 
yl=ql*sqrt (y+w0) ; 
y2=-ql*sqrt (y-w0) ; 
z1l=t-y2; 
z2=-yl; 
Z3=Z1/ ((z1*2)+(z2%2) ); 
z4=-22/ ((z1*2)+(z2%*2)); 
zl=t+y2; 
z2=yl; 
z5=(z1*z3) -(22*z4) ; 
z6=(z1*24)+(z2*z3); 
e(i)=z5; 
£(i)=z6; 
end; 
% 
% Find Pressure Contribution for Lower 
% 
21=0; 22=0;23=0; 24=0; z5=1; z6=0; 
for i=1:1:i1, 
zi=e(i); 
z2=£ (i); 
Z3=Z5; 
z4=z6; 
z5=(21*2Z3) -(z2*z4) ; 
26=(22*z3) +(z1*z4) ; 
end; 
z1=z5; 
22=26; 
t=t4*r9(n); 
z3=cos(t); 
z4=-sin(t); 
z5=(Z1*z3) -(z22*z4); 
z6=(22*23) +(Z1*z4); 
pl=pl1+ (s2*z5/r9(n)); 





p2=p2+ (s2*z6/r9(n)); 
% Truncation Test 
qa=qa* sqrt ((z5*z5)+(z6*z6)); 
if qa<0.00000001, break, end, 
end; 


% 

% Total Pressure at the Reciever 

% 

t5=sqrt (p1*2+p2%2) % Pressure at Reciever 
% Phase Angle at the Reciever 

phaseang=atan2 (p2,p1) *180/pi % Phase Angle 
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